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Abstract 

A recently proposed unified scaling law for interoccurrence times of earth- 
quakes is analyzed, both theoretically and with data from Southern California. 
We decompose the corresponding probability density into local-instantaneous 
distributions, which scale with the rate of earthquake occurrence. The fluctu- 
ations of the rate, characterizing the non-stationarity of the process, show a 
doubly power-law distribution and are fundamental to determine the overall 
behavior, described by a double power law as well. 
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Earthquakes constitute an extremely complex phenomenon in nature, with the deforma- 
tion and sudden rupture of some parts of the Earth crust driven by convective motion in the 
mantle, and the radiation of energy in the form of seismic waves. Only a part of this com- 
plexity is collected by earthquake catalogs, where magnitude, epicenter spatial coordinates, 
and starting time of events, among other measurements, are recorded. This information, 
which converts the phenomenon in a spatio-temporal point process marked by the magni- 
tude, nevertheless reveals some important scale-invariant properties. 

First, the Gutenberg- Richter law states that the number of earthquakes in some region 
with magnitude larger than some threshold value decreases exponentially with the threshold. 
Taking into account that (to a first approximation) the released energy increases exponen- 
tially with the magnitude, the probability distribution of the released energy turns out to 
be a power law, precisely the hallmark of scale-free behavior [1-3]. Second, the introduction 
of fractal geometry soon led to the recognition that the spatial distribution of epicenters (or 
hypocenters) draws a fractal object over the Earth surface [2,3]. And third, the Omori law, 
proposed more than 100 years ago, accounts for the number of events (called aftershocks) 
that follow a large shock after some time. This number is another power law, with exponent 
close to minus one [4] . 

This lack of characteristic scales suggests that the crust is in a critical state, like the well- 
known critical points studied in equilibrium systems, but without external adjustments of 
control parameters. Therefore, one may talk about a self-organized critical (SOC) state for 
the seismic system [5]. This concept has important implications for the issue of earthquake 
prediction [6]; indeed, a critical crust implies that a fracture process may or may not develop 
to provoke a large earthquake depending on minor microscopic details that are intrinsically 
out of control. 

An important quantity characterizing earthquake occurrence is the time interval between 
successive earthquakes. This time (that can be referred to as interoccurrence time, recurrence 
time, or waiting time) although related to the Omori law, has a distribution that is not 
clearly known. In fact, all the possibilities have been proposed, from periodic behavior for 
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large earthquakes to totally random occurrence. The most extended view is to consider the 
existence of two separated processes, one for the main shocks, which should occur randomly 
following a Poisson distribution, and another process for the aftershocks; but this should not 
hold for large events, for which clustering has been reported [7]. In general, the usual studies 
proceed by fixing a limited area of observation where aftershocks are skilfully identified 
and removed from data. On the opposite side, other works concentrate only on series of 
aftershocks. 

Bak, Christensen, et al. [8] have followed an alternative approach, which is to consider 
the problem in its complete spatio-temporal complexity. They divide the area of South 
California into regions of size L degrees in the north-south (meridian) direction and L degrees 
as well in the east-west (parallel) direction [9]. Only earthquakes with magnitude m larger 
than a threshold value m c are taken into account (but no other events are eliminated, all 
shocks are equally treated). For each L x L— region the time interval r between consecutive 
earthquakes is obtained for the period from 1984 to 2000 as t { = t» — i, where U is the time 
coordinate of the i—th earthquake within the region with m > m c . The probability density 
for this interoccurrence time, D(r,m c , L), is computed and the results give D(t, m c , L) oc 
1/r for short times and a faster decay for long times, with a dependence also on L and m c . 

Remarkably, when a scaling analysis is performed, all the distribution functions corre- 
sponding to different values of L and m c collapse into a single curve if the axis are rescaled 
by S b /L d f, with d f ~ 1.2, b ~ 1, and S = 10™ c (related to the energy roughly as S oc £ 2/3 ). 
In mathematical words, 

L d f ( L d f \ 1 ( L d f \ 

this scaling law constitutes the Bak-Christensen-Danon-Scanlon (BCDS) proposal [8]. For 
short times, the function G shows a slow variation not affecting the power-law (1/r) behavior; 
for long times, a fast decay is obtained, which could be consistent with an exponential 
distribution and therefore with a Poisson process, according to Ref. [8]. 

(From now on, to simplify the notation, we will omit the dependence of D on L and m c 
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and just write D{f).) 

This result is relevant for several reasons, among them: it shows scaling in the spatio- 
temporal occurrence of earthquakes, a key element to consider earthquakes as a critical 
phenomenon. Second, it is the first law that relates interoccurrence times, the Gutenberg- 
Richter law (factor l/S b ), and the fractal dimension of the spatial distribution of events 
(df), allowing a unified description. Third, the law is valid for all earthquakes, no matter 
their size or location, and no matter also if they are considered as aftershocks, foreshocks, 
or main shocks. Fourth, the power law tell us that immediately after any earthquake there 
is a high probability of having another one, and this probability decreases in time with no 
characteristic scale up to S b /L d f; that is, there is a correlation time that depends on the 
region size and magnitude under consideration, and therefore for any event one may find 
clusters of aftershocks in all time scales up to an appropriate length scale L. 

The importance of this law deserves further study. Here we are interested in a general 
understanding of the BCDS law and its origins. We will analyze the same catalogs as Bak 
et al. [10] and will show that the fast decay for long times is not exponential, but another 
power law. D{r) is related to its local and instant components and to the rate of earthquake 
occurrence r; this quantity, which counts the number of events per unit time in a given 
region, displays large fluctuations across several orders of magnitude, doubly power-law 
distributed. This is in contrast to simple SOC models. 

Let us pay more attention to the obtaining of the distribution D{r) by Bak et al. As 
we have mentioned, this distribution accounts for the time difference between successive 
earthquakes with magnitude larger than m c in every L x L— region. Times from different 
regions are counted together in D{r). But the total number of earthquakes differs from 
region to region (as it is well-know, due to the fractal spatial distribution), with a high 
variability [11]. Therefore, the local distributions D xy (r) accounting for the time difference 
in a given L x L— region (of spatial coordinates x, y) are clearly different. This means that 
D{t) is a mixed distribution constructed from all the different D xy {r). 

But further, looking into a single L x L— region one can see a high variability in the rate 
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along time t [11], see Fig. 1. In fact, the rate typically exhibits a quite stable behavior for 
some periods of time, with small fluctuations, but for other periods develops sudden burst 
of activity where its value increases sharply and then decreases to become stationary again, 
or not. This intermittency, of course related to the occurrence of larger earthquakes in the 
region, recalls the punctuated-equilibrium behavior of SOC systems [5,12], but note however 
that the variable that displays punctuated equilibrium is not only the signal m(t), but also 
the rate. 

Again then, the local distribution D xy {r) is obtained as a mixture of distributions, the 
densities of interoccurrence times in a given region at a certain time t, D xyt (r). From this, 
we can write 

D(t)kJ2 J D xyt (r)r(x,y,t)dt (2) 

\/x,y 

where the rate r(x,y,t) is the number of earthquakes per unit time in the region (x,y) at 
time t. (All functions here depend as well on L and m c , though the dependence is not 
explicitly written.) The rate r in the integral acts as a weight factor, due to the fact that 
the higher the rate in a given region and time, the larger the number of earthquakes that 
are produced and contribute to the distribution. 

We now make the hypothesis that the dependence on space and time enters into the 
distribution D xyt only through the rate r(x,y,t). That is, we assume that different regions 
at different times but with the same rate of occurrence will have the same distribution of 
interoccurrence times (if the rate is stationary), i.e., 

D xyt (T) = D(T\r(x,y,t)), (3) 

which is a conditional density. Therefore, 

D(t) = r D{r\r) r -^-dr, (4) 
J o /i 

with p(r) the probability density of the rate and \i = (r) just a normalization factor. 

Since D xyt (r) is an instantaneous quantity (and we have a single realization of the pro- 
cess), it were impossible to measure if we would not have the periods of stationarity in r. 
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Figure 2 shows these distributions for several periods of stationarity and several regions 
of different L and spatial coordinates. Indeed, the distributions D xyt do not only depend 
exclusively on r, but they scale with it, i.e., 

D(r\r) ~ r/(rr) (5) 

(a behavior that could have been derived by dimensional analysis) , with the scaling function 
/ being a power law for short times and having a fast decay for large ones. In fact, the 
distributions can be fitted by a function of the type 

f(u) = cJ^e-^\ (6) 

with 7 ~ 0.63, 5 ~ 0.92, u ~ 1.5, and C ~ 0.5. We could approximate then D(r\r) to a 
gamma distribution (5 = 1), which ensures that the large-scale cutoff is close to exponential. 

Next step is to look at the distribution of rates p(r). To be precise, r(x, y, t) is defined by 
counting the number of events above the threshold m c into the L x L— region of coordinates 
x, y during a time interval (t, t + At) and dividing the result by the duration of the interval, 
At. The corresponding probability density, which is calculated from 1984 to 2001 only for 
the regions in which there is earthquake activity, depends on L and At (as well as on m c ) and 
is shown in Fig. 3. In fact, the value of At should be small enough to ensure r ~constant, 
(but large enough for statistical significance), but there is no characteristic scale for constant 
r and therefore no typical value for At. Also, it is noteworthy in the figure the two-power-law 
behavior, one power law for low rates and another one for high rates, which can be modeled 
as 

P{r) = CO [ } a+p , (7) 
[1 + (6r) c }~ 

which gives p oc l/r 1_a for r <^ O^ 1 and p oc l/r 1+l3 for r ^> O^ 1 . We obtain exponents 
for r about 1 and 2.2, so a ~ and (3 ~ 1.2. Parameter c just controls the sharpness of 
the transition from one regime to the other and 6~ l is a scaling parameter. The fact of 
having two power laws means that there is no characteristic occurrence rate up to the value 
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O^ 1 , but for values in the tail of the distribution there is also scale invariance. This could 
be understood as criticality, not only in the time domain as we knew, but also in the rate 
domain. 

Additionally, it is easy to obtain the form of the scaling factor d^ 1 . The mean rate 
p = (r) is given by the total number of events divided by the total time and by the number 
of regions with activity; the former, because of the Gutenberg- Richter law scales as 1/S b 
and the latter as 1/ L d i , which gives (r) oc L d f / S b . Since the distribution turns out to scale 
in the same way, 0- 1 oc L d f/S b . For the scaling plot in Fig. 3 we have used b = 0.95 [8] 
and df = 1.6, which was obtained from a box-counting method for spatial distributions of 
epicenters with m > 2. 

Now that we know the form of the functions D(r\r) and p{r) we can answer the question 
about how the large variations of the rate influence the distribution D(t), just by integrating 
Eq. (4) with the use of (5)- (7). The limit r — > oo is obtained directly with the use of 
Laplace's method to evaluate asymptotic integrals [13]. We get 

D(t) oc — r'- 1 / r °+7 e -(Wu ) s dr oc — %— ; (8) 

P J0 5pT 2+a 

this is in fact independent on the tail of p(r), it does not matter if it is a power law or not. 
But there is another limit to be studied. Indeed, the behavior of the integral depends on the 
relation between r and 9. We have just calculated what happens for r 3> 9; the opposite 
case, r <^ 9, can be obtained for long times if we first perform the limit 9 — > oo and then 
apply Laplace's method for r. So, p{r) ~ C'9^ 13 /r 1+/3 , and the integral gives 

r 7-l roo r I 

Since p is proportional to 9~ l the last results can be summarized as follows, 

01-/3 01+a 

D(t) oc for r < 9, and D(r) oc for r > 9, (10) 

and are displayed in Fig. 4, where the two exponents for r turn out to be about 0.9 and 
2.2, giving (5 ~ 1.1 and a ~ 0.2 in good agreement with our calculations. Both results 
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for D(t) do not depend on the form given by Eq. (7) for p{r) as long as it exhibits the 
two-power- law behavior. In addition, the exponent 7 of D{r\r) does not affect the value of 
the two exponents of D{r). Notice that the scaling factor for r in D{r) is 9, that is, the 
inverse of the scaling factor for r, so, 9 oc S b /L d f. 

Finally, we would like to point out that D{r) is described by the same function (includ- 
ing the same values of the exponents) than the one that characterizes the trapping time 
distribution in a rice-pile model, see Fig. 1 of Ref. [14]. Also the coincidence between our 
Fig. 3 and Fig. 3 in Ref. [15] is notable, although no probability density is measured there. 
The exponents about 0.9 and 2.2 should be quite universal. 

In conclusion, we have performed a "microscopic" analysis of the BCDS law for earth- 
quakes, which provides a way to deal with the heterogeneity and non-stationarity of seismic 
occurrence. 

The author is profoundly indebted to Per Bak, who opened so many paths in Science, 
not only for his scientific guide, but for his personal warmth as well. Regarding this paper, 
he also thanks M. Boguha, K. Christensen, and Ramon y Cajal program. A fruitful part of 
this work was accomplished at l'Abadia de Burch (Pallars Subira, Lleida). 
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FIGURES 




1984 1986 1988 1990 1992 1994 1996 1998 2000 2002 
t (years A.D.) 

FIG. 1. Rate of earthquake occurrence as a function of time in the L = 10° region of South 
California, for m c = 2 with At = 2 months and for m c = 3 with At = 4 months. The vertical 
log-scale should not make underrate the large variations in r. Notice for example the constant rate 
at 1991 in contrast to 1992. 
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FIG. 2. Local distributions of interoccurrence times for several stationary periods and different 
regions, after scaling by the rate. The regions are labeled from x, y = to 10° /L — 1 from west to 
east (x) and from south to north (y). The fit is explained in the text; deviations at small times 
should be due to short scale disturbances of the stationarity. As in Bak et aVs paper, times smaller 
than 38 s are not considered. 
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FIG. 3. Scaled distributions of rates, for several At, L, and m c , using df = 1.6 and b = 0.95. 
Two power laws with exponents 1 and 2.2 fit the data. 



12 



0.01 



10" 



10" 



io- 



10" 



io- 



■ 'i \ ■ 'i ■ 1 t 1 ■ _r r 




L = 2.5°. rn c = 2 
L = 1.25°, ?n r = 2 
L = 0.625°, m c = 2 
i = 0.3125°, m c = 2 
L = 0.15625°, m c = 2 
L = 1.25°, m, = 3 
L = 0.625°, m r = 3 
L = 1.25° m„ = 4 
0.06/m - 9 
7000/m 2 - 2 

J lJ I I J I 



0.01 0.1 



10 100 1000 10 4 



10 : ' 



10° 



10' 



FIG. 4. Scaled distributions of interoccurrence times, D(t), for different L and m c with dt = 1.6 
and b = 0.95; r > 38s again. The straight lines illustrate the double power-law behavior, with 
exponents 0.9 and 2.2. 
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